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Abstract 

Recent advances in optical systems make them ideal for undersampling 
multiband signals that have high bandwidths. In this paper we propose 
a new scheme for reconstructing multiband sparse signals using a small 
number of sampling channels. The scheme, which we call synchronous 
multirate sampling (SMRS), entails gathering samples synchronously at 
few different rates whose sum is significantly lower than the Nyquist sam- 
pling rate. The signals are reconstructed by solving a system of linear 
equations. We have demonstrated an accurate and robust reconstruction 
of signals using a small number of sampling channels that operate at rel- 
atively high rates. Sampling at higher rates increases the signal to noise 
ratio in samples. The SMRS scheme enables a significant reduction in the 
number of channels required when the sampling rate increases. We have 
demonstrated, using only three sampling channels, an accurate sampling 
and reconstruction of 4 real signals (8 bands). The matrices that are used 
to reconstruct the signals in the SMRS scheme also have low condition 
numbers. This indicates that the SMRS scheme is robust to noise in sig- 
nals. The success of the SMRS scheme relies on the assumption that the 
sampled signals are sparse. As a result most of the sampled spectrum 
may be unaliased in at least one of the sampling channels. This is in 
contrast to multicoset sampling schemes in which an alias in one channel 
is equivalent to an alias in all channels. We have demonstrated that the 
SMRS scheme obtains similar performance using 3 sampling channels and 
a total sampling rate 8 times the Landau rate to an implementation of 
a multicoset sampling scheme that uses 6 sampling channels with a total 
sampling rate of 13 times the Landau rate. 

1 Introduction 

In many applications of radars and communications systems it is desirable to 
reconstruct a multiband sparse signal from its samples. When the carrier fre- 
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qucncies of the signal bands are high compared to the overall signal width, it 
is not cost effective and often it is not feasible to sample at the Nyquist rate. 
It is therefore desirable to reconstruct the signal from samples taken at rates 
lower than the Nyquist rate. There is a vast literature on reconstructing signals 
from undersampled data [5] — [5]. Most of methods are based on a multicoset 
sampling scheme. 

In a multicoset sampling scheme m low-rate cosets are chosen out of L cosets 
of samples obtained from time uniformly distributed samples taken at a rate F 
where F is greater or equal to the Nyquist rate F nyq [3J. In each channel the 
sampling is offset by a different predetermined integer multiple of the reciprocal 
of the rate F. The data from the different sampling channels are then used 
to reconstruct a signal by solving a system of linear equations. Under certain 
conditions on the sampling rate and number of channels, a proper choice of the 
time offsets ensures that the equations have a unique solution in case that the 
signal bands are known a priori [3J , or unknown a priori [5] . 

A previous work has demonstrated a different scheme for reconstructing 
sparse multiband signals [9]. The scheme, called multirate sampling (MRS), 
entails gathering samples at P different rates. The number P can be small 
and does not depend on any characteristics of a signal. The approach is not 
intended to obtain the minimum sampling rate. Rather, it is intended to re- 
construct signals accurately with a very high probability at an overall sampling 
rate significantly lower than the Nyquist rate under the constraint of a small 
number of channels. 

The reconstruction method in 9] does not require synchronization of the 
different sampling channels. In addition to reducing the complexity of the sam- 
pling hardware, unsynchronized sampling relaxes the stringent requirement in 
multicoset sampling schemes of a very small timing jitter in the sampling time 
of the channels. Simulations have indicated that MRS reconstruction is robust 
both to different signal types and to relatively high noise. 

Accurate signal reconstruction using the unsynchronized MRS scheme of [5] 
requires that each frequency in the support of the signal be unaliased in at least 
one of the sampling channels. In this paper we describe a new reconstruction 
algorithm that overcomes this deficiency by using synchronized sampling chan- 
nels. In our synchronized multirate sampling (SMRS) scheme, the aliasing is 
resolved, in the same spirit as multicoset sampling schemes, by solving a system 
of linear equations. Therefore, the SMRS scheme enables the reconstruction of 
signals that cannot be reconstructed using MRS scheme. 

In our SMRS scheme each sampling rate must be an integer multiple of 
the same basic frequency. The reconstruction method in our SMRS scheme 
requires the same frequency resolution in all sampling channels despite the fact 
that the sampling rate is different in each channel. This requirement is well 
suited to the implementation of the scheme using the optical system described 
in [TU]. The sampling is performed in two steps. In the first step the entire 
signal spectrum is downconverted into a low frequency region called baseband 
by multiplying the signal with a train of short optical pulses [10]. In each of 
the sampling channels the pulse rate is different. The frequency downconverted 
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analog signals are then converted in each sampling channel to digital signals 
using an A/D electronic converter that samples at the highest of the channel 
rates. The result is that the signal is sampled in all channels with the same 
time resolution that is determined by the sampling rate of the A/D converter. 
Alternatively, since each sampling rate must be an integer multiple of the same 
basic frequency, a common frequency resolution can be obtained by using the 
same time window for all channels. 

There is an inherent advantage to sampling, in each channel, near the max- 
imum sampling rate made possible by cost and technology. This is because 
sampling at higher rates increases the signal to noise ratio in sampled signals 
[9 . Our simulation results indicate that when the sampling rate in each channel 
increases, our SMRS scheme requires significantly fewer sampling channels than 
does a multicoset sampling scheme of [5] to obtain comparable reconstruction 
success. When the sampling rate in each channel increases, the probability that 
a sparse signal aliases simultaneously in all sampling channels becomes very low 
in an MRS scheme. It is lower than in a multicoset sampling scheme in which, 
because all channels sample at the same frequency, an alias in one channel is 
equivalent to an alias in all channels. Our numerical simulations indicate that 
the success rate of our SMRS scheme is significantly higher than that of a mul- 
ticoset sampling schemes of [3J when the number of sampling channels is small 
(3 in our simulations) and the sampling rate of of each channel is high. 

Exactly the same data obtained in the SMRS scheme can be obtained using 
a multicoset sampling scheme. However the multicoset pattern requires many 
more channels, each of which samples at a very low rate. In a numerical example 
of section l4~T1 a 3 channel multirate sampling pattern is equivalent to a 58 channel 
multicoset sampling pattern. With this multicoset pattern the time difference 
between two consecutive samples is on the order of 1 psec. Such an accuracy 
cannot be practically achieved. 

In an SMRS scheme the data is reconstructed differently than in a multicoset 
scheme implementation of [SJ . In this multicoset recovery scheme it is assumed 
that, in addition to being sparse, the spectrum of the signal consists of bands 
each of which is narrower than the coset sampling rate. However, our SMRS 
scheme requires no such assumptions on the originating signal. The reconstruc- 
tion of sparse signals in our SMRS scheme is robust. This is because most of the 
sampled spectrum is unaliased in at least one of the sampling channels. This 
is in contrast to the equivalent multicoset sampling scheme with many low rate 
sampling channels in which the alias probability is very high and an alias in one 
channel is equivalent to an alias in all channels. In such cases a blind signal 
recovery of [5J can be found only by using a pursuit algorithm whereas, in many 
cases, the SMRS scheme doesn't require a pursuit algorithm. The signal is re- 
constructed directly by a single matrix inversion. By making a very reasonable 
physical assumption, we are also able to simplify the reconstruction by reducing 
the number of possible signal locations in a straightforward manner. 

The SMRS scheme reconstructs signals by solving a system of linear equa- 
tions. Our simulations indicate that linear equations in an SMRS scheme are 
numerically stable in that they have low condition numbers. This makes them 
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rather insensitive to noise. In the SMRS scheme, when sampling at total sam- 
pling rate that is significantly higher than theoretically required, the probability 
that part of the signal spectrum will not be aliased in at least one of the sam- 
pling channels is very high. The unaliased parts of the signal can be recovered 
directly from the sampling channels. Our simulation results indicate that the 
sensitivity of the reconstructed signal to noise added at frequencies of baseband 
that are not aliased in at least one of sampling channels is very low. 



2 Synchronized MRS 

In this section we describe the SMRS scheme. Let F max be an assumed maxi- 
mum carrier frequency and let X(f) S L 2 ([0, F max \) be the Fourier transform 
of a complex-valued signal x(t) that is to be reconstructed from its samples. 
Throughout the analysis we normalize the Fourier transform by 



X(f) = 




dt. 



The modifications required to reconstruct real-valued signals are described in 
the Appendix. We assume that the signal x(t) to be sampled, in addition to 
being bandlimited in [0, F max ], is multiband; i.e., the support of its spectrum 
is contained within a finite disjoint union of intervals (a n ,b n ], each of which 
is contained in [0, i^ max ]. By assumption, max6„ < F max . In reconstructing a 
signal we do not assume any a priori knowledge of the number or location of the 
intervals (a n ,b n ]. The only requirement is that the signal be sparse; i.e., that 
for a signal whose spectral support is contained within the N intervals (a, , bi] , 

J2k=i b k ~ a k < -Fmax- 

In the MRS scheme of [5] a signal is sampled at P different sampling rates 
F % (i = 1 . . . P) . If the delays of the channels arc denoted by A 1 , the sampled 
signals x l (t) are given by 

oo 

x i (t)=x(t) <K*-|?- Ai ) « 

71— — OO 

where S(t) is a dirac delta "function". The spectrum X l (f) of the sampled 
signal in the ith channel satisfies 

oo 

X\f) = F i X(f + nF*)exp{j2Tr(f + nF*)A*}. (2) 

71 — — OO 

Since the channels are synchronized in time, we may set each A* = 0. Equation 
([21) becomes 

oo 

X i (f)=F i X{f + nF l ). (3) 
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It follows from ^ that all the information about the ith sampled spec- 
trum X l (f) is contained in the interval [0, F 1 ]. This interval is called the ith 
baseband. To process the sampled signals it is necessary to discretize the fre- 
quency axis. We use the same frequency resolution (Af) in all of the sampling 
channels. The same frequency resolution is directly obtained if the system is 
implemented using an optical sampling system [lOj . To use our reconstruction 
method each sampling rate should be chosen to be an integer multiple M 1 of 
the basic frequency resolution: F l = APAf. By using a sampling time window 
with a duration T = l/Af, the same frequency resolution Af is obtained in 
all of the sampling channels. Alternatively, the same frequency resolution can 
be obtained in an optical system by using A/D converters with the same rate 
at the output of all the optical downconvertors [10] • To reduce computational 
requirements we ignore redundant data. Thus, only AP entries of the Discrete 
Fourier transform (DFT) are retained for further processing. 

We represent the signals over the common discretized frequency axis with 
the following notations: 

X l [k] = X l (kAf), fc = 0,...,M l -l, 
X[k] = X OA/), fc = 0, ...,M-1, 

where Af is the frequency resolution, M = \F max /Af~\, X l [k] is the spectrum 
of the sampled data from the discretized frequencies in the baseband [0, F l ] and 
X[k] is the spectrum of the originating signal at the discretized frequencies. 
Equation ([3]) takes the form 

M-l oo 

X i [k] = F i ^2 5[l-{k + nM 1 )] (4) 

£=0 n=— oo 

where S[n] is the Kronecker delta function. Equation Q can be written in 
matrix form as follows: 

x J = A l x (5) 

where x l and x are given by 

(x^fc+i =X i [k], 0<fc<ikP-l, (6) 
(x) fe+1 =X[k] 0<k< M-l, 

and A 1 is a M % x M matrix whose elements are given by 

oo 

K+u+i=F l S[l-(k + nM% (7) 

n=— oo 

Each element A l k l is equal to either F l or 0. This is because there is at most 
one contribution in the infinite sum of <5's which is made when I = k( mod M l ). 
In each row of the matrix A 1 there are only [F max /F n \ non zero elements. 

For each value of i (i = 1 . . . P) , ([5]) defines a set of linear equations that 
relate the spectrum of the signal to the spectrum of its samples. The vector 
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x in iJSJ) is the same for all the P equations because it doesn't depend on the 
sampling. 

The vector x and the matrix A are obtained by concatenating the vectors 
x* and matrices A 1 as follows: 



/ x 1 \ 




f A l \ 


x 2 




A 2 




, A = 











These form the system of equations 

x = A x. (8) 

The matrix A has exactly P non-vanishing elements in each column that cor- 
respond to the locations of the spectral replica in each channel baseband. 
In case that the signal is real-valued its spectrum fulfills 

X(f)=X(-f) (9) 

where a + bj — a—bj is the complex conjugate. The equations for reconstructing 
such a signal are described in the Appendix. 

To invert ((8]) and calculate the discretized signal spectrum (x) it is necessary 
that the number of rows Yl^—i M l in A be equal to or larger than the number 
of columns M. Defining Ftotal = Sj=i F l makes this condition equivalent to 
the condition 

-Ftotal > fmax- (10) 

The condition on the sampling rates given in p0[) is consistent with the require- 
ment that the sampling rate be greater than the Nyquist rate of a general signal 
whose spectral support is [0,F max ]. However, when sampling sparse signals, an 
inversion of the matrix may be possible even if the condition (|10[) is not fulfilled. 
Our objective is to invert JHJ) in the case of sparse signals with sampling rates 

-^total ^ -Fmax- 

3 Inversion Algorithm 

In this section we describe our inversion algorithm for the SMRS scheme. The 
purpose of the algorithm is to invert (jHJ; i.e., to calculate the vector x from 
the vector x. To invert the equations with sampling rates lower than those 
prescribed by (|10[) the assumption that the signal is sparse should be taken into 
account. 

3.1 Known Band Locations 

In the case in which the signal band locations (a n , b n ] are known ([8]) can be 
simplified easily. All the elements of x that correspond to the frequencies not in 
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the spectral support U„(a n , b n ] are eliminated from ©. All the columns of the 
matrix which correspond to these elements arc also eliminated. The reduced 
system of equations that corresponds to ([5]) is given by 

X r ed = -A. re dX rc d. (11) 

A unique solution exists only if A ro d is full column rank. In this case the 
inverse can be found using the Moore-Penrose pseudo- inverse [12] , In a matrix 
of full column rank the number of rows must equal or exceed the number of 
columns. 

Although the entire spectrum is downconverted to baseband, we assume that 
we are sampling highly sparse signals. Hence the number of non-zero entries in 
each x l is significantly smaller than M l in at least one of the sampling channels. 
The number of non-zero entries in each x 1 might even be smaller because of 
aliasing. A necessary condition for a unique inverse or pseudo-inverse of A rc d is 
that this number still be greater than the number of non-zero entries of x. This 
is consistent with a Landau theorem p] that states that one cannot reconstruct a 
signal if the spectral density of the samples collected from all sampling channels 
is less than the spectral support of the originating signal. 

The choice of sampling rates imposes restrictions on the possible values of 
-Fmax for which an inversion of (jlip is possible. For the matrix A rc d to have full 
column rank, it must not have any identical columns. Since we do not restrict 
the possible locations of the known signal bands, any combination of columns 
of the matrix A may appear in the matrix A ro( j. Therefore we require that A 
not have any identical columns. The matrix A is composed of P sub-matrices 
A' whose columns are periodic: 

For the matrix A not to be periodic, it is required that any common period 
of the P sub-matrices be larger than M . This condition is met if the least 
common multiple of the {M l }i is larger than M. As a result, F meix should fulfill 
-Fmax < lcm(M 1 , . . . , M p )Af, where lem denotes least common multiple. 

3.2 Unknown Bands' Location 

In the case that the locations of the bands (a n , b n ] are not known a priori some 
additional assumptions must be made. In the multicoset recovery scheme of [5] 
it was assumed that the maximum band width is known a priori. In our SMRS 
scheme we do not make any assumptions on the intervals (a„,6„] but instead 
we add assumptions on a signal's spectrum itself. 

We assume that, for each discretized frequency kAf, any sampled spectrum 
X l (kAf) — is due only to lack of a signal in any of its replicas and not due 
to any aliasing. In other words if, for any n, X(f + nF 1 ) ^ 0, then X l (f) ^ 0. 
Another assumption is that there is a unique solution in the case the signal 
support is known, i.e., the matrix A re d has a full column rank for a known 
support. 
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Applying the first assumption, one can detect baseband frequencies in which 
there is no signal. These frequencies can be eliminated in the reduction of ([8]). 
We describe this simple procedure for eliminating frequencies which, according 
to our assumption, cannot be part of the spectral support of the originating 
signal. The elimination is similar to one presented in asynchronous-MRS !) . 
We denote the indicator function x l [i] as follows: 

<m _ J 1, for all I G [0, M - 1] such that X l (lAf) ^ , « 

X [ 1 ~ \ 0, otherwise. { ' 

The function X l (f) is periodic with period F % . Therefore % l [Z] ^ s a periodic 
extension of an indicator function over the baseband / e [0, F l ). 
We define the x[l] as follows: 

p 

x[l] = l[x% ie[0,M-l]. (13) 

i=l 

The function %[/] equals 1 over the intersection of all the upconverted bands 
of the P sampled signals and it defines the columns of the matrix A that are 
retained in forming the reduced matrix A ro( j. All other columns are eliminated 
and their corresponding elements in the vectors x are also eliminated. After the 
elimination of the columns from the matrix A, the matrix rows which correspond 
to zero elements in x and their corresponding elements in the vectors x are also 
eliminated. In some cases the function \\f\ equals 1 only for frequencies within 
the spectral support of the signal. In such cases the resulting equations are 
identical to those found in the previous subsection (equation (fTTj) ). However, 
in other cases, x[Z] may also equal 1 for frequencies outside the signal's true 
spectral support. In such cases the reduced matrix will have more columns than 
the matrix in the case in which the spectral support of the signal is known. As 
a result the inversion requires finding the values of more variables. 

Each eliminated zero energy baseband component causes elimination of re- 
spective rows and columns. The elimination of one baseband entry means that 
all the frequencies that are downconverted to that baseband entry (the aliasing 
frequencies) are also eliminated. This is because of our first assumption that 
zero entry in the baseband corresponds to zero entries in all of the frequency 
components of the original signal that are down-converted to frequency of the 
baseband entry. Therefore, elimination of one baseband entry results in elimi- 
nation of Lfmax/ min{i ;i4 }J to \F max / max{F l }~\ corresponding columns. Thus, 
if the number of the zero elements in x is sufficiently large, the number of rows 
in the matrix A rec j may be larger than the number of columns. 

If in addition, matrix A rec j has a full column rank, the problem is cither 
consistent or overdetermined. In such cases there is a unique inversion to 
which can be found using the Moore-Penrose pseudo-inverse. If the matrix is 
not full column rank, the problem is underdetermined and the inversion is not 
unique. A unique solution in such cases can be obtained either by increasing 
the total sampling rate or by adding additional assumptions on the signal. 
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3.3 Ill-posed cases 



In many cases the matrix in for unknown band locations is not full column 
rank. In these cases there are subsets of columns in the matrix A re( j that are 
linearly dependent. Using this linear dependence, a solution to (TTTj) can be 
found. However any solution found can be used to construct an infinite number 
of solutions to the equation. Thus, there is no unique inversse to and the 
inversion problem is ill-posed. 

To reconstruct a signal in the case in which the inversion problem is ill- 
posed we impose an additional assumption on the signal. We assume that in 
the case the signal support is unknown and the problem is ill-posed, among all 
possible solutions to (fTT|) . the originating signal is the one that is composed of 
the minimum number of bands. This is the signal we attempt to construct. We 
also assumed earlier that the reduced matrix which corresponds to the case in 
which the signal support is known is well-posed. In case that this assumption is 
not fulfilled the sampled signal does not contain enough information for solving 
the problem. 

Under the three assumptions stated above (the assumption that leads to 
matrix reduction, the existence of the unique solution to (jlip when the signal 
bands are known, and band-sparsity) the inversion problem is reduced to finding 
the solution of (fTTj) that is composed of the minimum number of bands. The 
problem is NP-complex since we need to test every possible combination of 
bands. 

The algorithm described here is of lower complexity and its purpose is to find 
a solution of (fTTj) that is composed of the minimum number of bands without 
testing all the combinations. The resulting algorithm attains a lower success 
rate but decreases the run-time significantly as compared to an NP-complcx 
algorithm. We do not provide the conditions under which the correct solution 
is obtained. 

Our algorithm is based on the Orthogonal Matching Pursuit (OMP) [8]. 
This algorithm belongs to the category of the "Greedy Search" algorithms. The 
original OMP algorithm is used to find the sparsest solution x of underdeter- 
mined equations Ax = b [8] where A is an underdetermined matrix. The 
sparsest solution is the solution having the smallest norm ||x||o where ||x||o is 
the number of non zero elements in the vector x. The original OMP algorithm 
collects columns of the matrix A iteratively to construct a reduced matrix A r . 
At each iteration n the column of A which is added to A™ -1 to produce a matrix 
A™ is the column which results in the smallest residual error min x ||b — A"x||2 
where for every vector y, ||y||| = J2iVi- The iterations are stopped when some 
threshold e is achieved. Sufficient conditions are given for the algorithm to 
obtain the correct solution [8]. 

We denote A = A re d, b = x rc d and x = x ro d. Since we are seeking the 
solution of Ax = b with the smallest number of bands and not the smallest norm 
jx|| , we modify the OMP algorithm by instead of choosing a single column as 
in [5] by selecting iteratively blocks of possible locations. The columns of the 
matrix A fall into J blocks. Each block is indexed by j and Bj contains the 
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index of columns of the jth block. 

Each index set Bj identifies a possible band of the spectral support of the 
reconstructed signal. We start the iteration with the empty set S° = of 
column indexes, the empty matrix Aj?, and the set B° = Uj=i Bji so that at 
nth iteration the following holds: S n [J B n = B°. At the nth iteration (n > 1) 
the algorithm must decide which block to add to A™ -1 . If the index set Bj 
is chosen, then S n = S" 1 " 1 U Bj and B n = B^ 1 \ B r The matrix A? is the 
matrix whose columns are selected from A according to the indexing set S n . 

The block added is the one that produces the smallest residual error e™ = 
min Bj 6 ju-i min x ||b — A"jX||| where A"j is the matrix obtained by adding the 
block indexed by Bj to A" -1 . The algorithm stops when the threshold e is 
reached. The threshold e is a very small number and reflects upon the finite 
numerical precision of the calculations. 

The algorithm performed well in our simulations. However, there were cases 
in which the support of the reconstructed solution did not contain all the orig- 
inating bands and cases in which the reconstructed signal was incorrect even 
though all the assumptions on a signals given in section f3. 21 were fulfilled. Only 
after performing the last step of the algorithm is it possible to determine the 
support of the signal. The algorithm failed primarily for one of two reasons. 
One of them was due to the inclusion of a block that reduced the residual error 
on one hand but on the other hand, caused a resulting matrix A™ to be not 
full column rank as hypothesized in our problem (in section [JO]). This can hap- 
pen, for example, when a block consists of a correct sub-block and erroneous 
sub-blocks. Including any erroneous sub-blocks may result in an ill-posed prob- 
lem. Another reason for failure was a large dynamic range of the signals. When 
reconstructing such signals, correct bands may be ignored by the algorithm in 
cases that the energy within the bands is significantly lower than the energy in 
other bands. 

Sufficient conditions that assure that the algorithm converges to a unique 
solution have not yet been determined. 

4 Simulations results 

The ability of the signal reconstruction algorithm to recover different types of 
signals was tested. In one set of simulations the ability of the algorithm to 
reconstruct multiband complex and real-valued signals with different spectral 
supports, shapes, and band widths that were not known a priori was tested. 
Additional simulations were performed in which real-valued multiband signals 
were contaminated by additive white noise. Band carrier frequencies were chosen 
from a uniform distribution over the maximum support: 0-20 GHz for complex 
signal and —20 to 20 GHz for real- valued signals. For each set of simulations we 
counted the mean rate of ill-posed cases in which the modified OMP algorithm 
was used to recover the signal. Mean times for accurate signal reconstruction 
were also recorded. Failures of the reconstruction were either because one of 
the initial assumptions was not fulfilled or because of the failure of the modified 
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OMP algorithm. 

In different simulations the width of each band (hence the total bandwidth 
of the signal) was varied. The number of bands was always set equal to 4 for 
complex signals and to 8 for real-valued signals (4 positive bands and 4 negative 
bands). However, the reconstruction algorithm was unaware of this number. In 
all the simulations the frequency resolution was set to 5 MHz. The simulations 
were performed on a 2 GHz Core2Duo CPU with 2 GB RAM storage in the 
MATLAB 7.0 environment (no special programming was performed to use both 
cores) . 

4.1 Ideal multiband signals 

In the case of ideal multiband signals, because of the absence of energy outside 
of strictly defined bands, one can expect a perfect reconstruction. Accordingly, 
the algorithm was evaluated by a perfect reconstruction criterion; i.e., a mean 
difference between the true and the reconstructed signal spectrum less than 
1CP 10 . Whenever this error was attained, the reconstruction was deemed to 
have been successful. Otherwise, it was deemed to have failed. The threshold 
for the modified OMP was chosen accordingly; e = 1CP 20 . 

Simulations were carried out for complex signals to compare the results to 
those published using the multicoset sampling recovery scheme of [5.. The 
sampling rates were chosen to be 0.95, 1.0 and 1.05 GHz yielding a total sampling 
rate F totai = 3.0 GHz. 

Different signals with 4 bands of equal width were generated. Each band 
was chosen to lie within the interval [0, 20] GHz. Both the real and imaginary 
spectra of the signal within each band were chosen to be normally distributed. 
Specifically, for each frequency / = kAf in a chosen band, the real and imag- 
inary components of X(f) were chosen randomly and independently from a 
standard normal distribution. Each bands' spectra were scaled by a constant 
a such that each bands' energy was equal to a uniformly generated value E on 
the interval [1,5]; i.e., for specific band, 

X(f) = a[X r (f) +jX im (f)], \\X(f)\\ 2 = E. 

These signals were also used to test the multicoset sampling reconstruction 
scheme of [5|. The empirical success rates were obtained from 1000 runs, each 
with a different total bandwidth (-FLandau)- The success rate is shown in Fig. [TJ 
We have validated that the empirical success rate did not significantly change 
when the number of simulation runs was increased from 1000 to 5000. 

As is evident from Fig. [TJ the empirical success percentage of an ideal recon- 
struction is high when Ftotai/FLandau > 5. In the SBR2 scheme (downsampling 
factor L = 199) in [5] the empirical success rate shows perfect reconstruction 
for at least p — 14 channels. This corresponds to F^tai = p/LT — 1.4 GHz and 
hence F^tai /FLandau should be greater than about 3. Although the total sam- 
pling rate in [5] is lower than in SMRS, the number of channels that are used in 
that scheme is significantly higher compared to that used in SMRS where only 
3 channels are used. 
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Another simulation in [5] shows that for lower number of channels with 
L = 23 empirical perfect reconstruction is achieved with at least 6 channels and 
-Ftotai/^Landau > 13. In the SMRS scheme empirical perfect reconstruction was 
obtained using only three channels with a total sampling rate -Ftotai/-FLandau > 5. 
The system parameters (number of sampling channels, sampling rates, F mlLX ) 
that were used in our last simulation are the same as those used in our optical 
sampling experimental set up. The fact that the simulation results were obtained 
in a practical situation demonstrates that our SMRS scheme can reconstruct 
sparse signals perfectly using both a fewer number of sampling channels and with 
a lower total sampling rate than are required by multicoset sampling schemes. 

We note that the same data that is obtained in a SMRS pattern can always 
be obtained by a multicoset sampling pattern since the ratio between each pair 
of sampling rates is rational. In our example the sampling rate of each coset 
is equal to 1/LT = 50 MHz. The number of multicoset sampling channels (p) 
is 58. The time offset between the cosets is a multiple of T — 399 ]jhz ■ ^hc 
downsampling factor L is 399 GHz/50 MHz=7980. Note that since L is not 
prime, we are not guaranteed to obtain a universal sampling pattern [5]. 

Fig. [1] shows that in our scheme 100 percent empirical success was obtained 
for -Ftotai/^Landau > 5. This corresponds to a maximum overlap in a coset 
scheme of K = 16 and hence p/K = 3.6 for p = 58. In comparison, for 
p = 14 channels and downsampling factor of L = 199 in a multicoset sampling 
scheme in [5], the empirical success rate of 100 percent was obtained for SBR2 
algorithm for p/N = 3.5. This ratio is similar to the equivalent ratio obtained in 
our scheme (3.6) although the downsampling factor in our equivalent multicoset 
scheme is equal to 7980. For approximately 95 percent success rate in SBR2 
scheme of [5], the ratio p/N is 2.75. In our scheme, for the same empirical 
success rate of 95 percent, F to t a i/-FLandau = 3.53 and the equivalent ratio of the 
number of channels to maximal number of overlaps is p/K = 2.4. 

The mean percentage of ill-posed cases is shown in Fig. [2j The figure shows 
that for 4 < Ftotai/^Landau < 10, in most of the tested cases the inversion was 
ill posed. Nonetheless, a very high success percentage was obtained for these 
values. This indicates that our modified OMP algorithm was very successful in 
resolving these cases. 

Fig. [3] shows the mean run time as a function of ^totai/^Landau (constant 
total sampling rate and varying signal support). Because matrix inversion is 
the most computationally intensive operation in the algorithm, mean run time 
decreases with a reduction in signal bandwidth. This is because, with a fixed 
resolution, matrix size is proportional to signal bandwidth. Moreover, as the 
ratio -F to tai/-FLandau increases, the possible spectral support obtained at the first 
step of the reconstruction increases beyond the increase of the signal bandwidth. 

We note that we could reduce the run time without significantly affecting 
the empirical success percentage by solving ((8|) using a low resolution and re- 
constructing the signal using a higher resolution. 

The algorithm, modified as explained in the Appendix, was also tested 
against real- valued signals. The assumed maximum frequency F max was set 
to 20 GHz. The number of sampling channels was set to 3 with the sampling 
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Figure 1 : Empirical success percentages for 4-bands complex signals for different 
spectral supports (.FLandau) with i*kyquist=20 GHz, total sampling rate -Ftotai = 
3 GHz 
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Figure 2: Ill-posed cases mean percentage for 4-bands complex signals for dif- 
ferent spectral supports (-FLandau) with -FNyquist=20 GHz, total sampling rate 
-Ftotai = 3 GHz 



frequencies chosen to be F 1 = 3.8 GHz, F 2 = 4.0 GHz, and F 3 = 4.2 GHz re- 
sulting in a total sampling rate, Ptotal — 12 GHz. The sampling frequencies are 
the same as are used in our experiments based on asynchronous MRS [SJ. The 
number of bands was set to 8 (4 positive and 4 negative frequencies, assuming 
no carrier frequency so low as to have the frequency in its spectrum). Each 
band was chosen to be of equal width -FLandau/8. 

Once a band (a, b] was chosen, the spectrum of X(f) for / 6 (a, b] was 
determined by the following formula: 



X(f)=Asm 



7r(/-a) 
b — a 



The phase 9 was chosen randomly from a uniform distribution on [0, 27r] and 
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Figure 3: Mean run times for 4-bands complex signals for different spectral 
supports (^Landau) with i<Nyquist=20 GHz, total sampling rate F to tai = 3 GHz 



the amplitude A was chosen randomly from a uniform distribution on [1, 1.2]. 

Fig.0]shows the empirical success percentages of the algorithm tested against 
real valued signals. As is evident from the figure, the empirical success percent- 
age is high when -Ftotai/-fLandau > 8. We note that the required sampling rate 
is significantly higher in this example than in the complex signals simulation. 
The reason is that in our real case example there are twice as many bands as in 
the complex case simulations. Hence, after the sampling, an overlap may also 
occur between the negative and the positive bands of the real signal. We note 
that when sampling a real signal at a sampling rate F l , it is sufficient to know 
the spectrum in a frequency region [0, F l /2}. However, for real signals, there is 
uncertainty as to whether a signal in baseband is obtained from a signal in the 
positive band or in the negative band. 

The number of ill-posed cases and the mean recovery run times for the real- 
valued signals are shown in Figs. [5] and [6] respectively. It can bee seen that the 
mean rate of ill-conditioned cases is much lower for real- valued signal simulations 
than for complex ones. This could be due to the correlation between positive 
and negative frequency components of real signals. 

4.2 Solution Stability 

The stability of the linear equations used in the recovery scheme was tested 
via the condition number for the real- valued signals equations f (|2ip in the Ap- 
pendix). This case is important in our experiments since we sample real sig- 
nals [9] . The reconstruction scheme for real- valued signals requires solving two 
systems of linear equations; one for the real part and the other for the imaginary 
part. Each of the two systems of equations is described by a different matrix. 
In each test case we presented the maximum value of the condition numbers 
of the two matrices. For 4-band 200-MHz-width randomly generated signals 
(-FLandau=L6 GHz) the condition number among 1000 runs was at most 5.3. 
The condition numbers histogram is shown in Fig. [Jj 
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Figure 4: Empirical success percentages for equal 4-bands real signals for dif- 
ferent spectral supports (-FLandau) with i*Nyqmst=40 GHz, total sampling rate 
Ftotal - 12 GHz 
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Figure 5: Ill-posed cases mean rate for 4-bands real signals for different spectral 
supports (F Landau ) with F Nyquist =40 GHz, total sampling rate F tot al = 12 GHz 



When sampling a sparse signal, most frequencies of the signal are unaliased 
in at least one of the sampling channels. It can be easily shown that when a 
frequency component of the signal is not aliased in any of the sampling chan- 
nels, the reconstructed signal is obtained simply by averaging the corresponding 
sampled spectrum at the different sampling channels. The reconstruction of a 
frequency component that is unaliased in at least one of the sampling channels 
can also be easily performed by copying the corresponding unaliased spectrum 
to the reconstructed signal. Therefore, for sparse signals, the reconstruction in 
the SMRS scheme is robust and the condition number of the matrices is small. 
Indeed, we have verified that the condition number of the matrix increases as the 
number of frequencies in the original signal that are aliased increases. Fig. [8] 
shows the mean values of condition numbers versus number of aliased frequen- 
cies as calculated for the 4 real signals each with a 200 MHz width that are 



15 



sampled at F to tai = 12 GHz. The other parameters of the simulation are the 
same as those in the simulation that resulted in Fig. |4j 

4.3 Noisy signals 

The algorithm's performance was also tested for its ability to reconstruct real- 
valued signals contaminated by Gaussian white noise. The presence of noise 
demands some modification of the algorithm. One modification is in detecting 
the possible bands of the support of the originating signal. Because the spectral 
support of white noise is not restricted to the spectral support of the uncontam- 
inated signal, the indicator functions in (|12p cannot be used. Instead, we adapt 
(jX3J) to noisy cases similarly as in [5]. In [S], for the indicator function x l [l] to 
be equal to 1 at any frequency, it was required that the average energy of the 
signal in the neighborhood of that frequency be higher than a certain threshold. 
In SMRS we further expand each band in x[l] to include additional frequencies 
that might otherwise be omitted when defining the indicator functions X*[/]. 
Once the bands are identified the matrix equations are constructed exactly as 
in the noiseless case. 

The solution of the linear equations given in (|21[) is modified in the noisy case. 
Because the added white noise affects the entire spectrum, a signal contaminated 
by white noise can no longer be considered multiband in the strict sense. Thus 
one cannot expect to reconstruct it perfectly from samples taken at a total 
rate lower than the Nyquist rate. Whereas in the ideal noiseless case the error 
norm vanishes, with a signal containing noise, one must relent on a perfect 
reconstruction and settle for a minimum error. In the noisy case the solution 

to (f2"Tj) should solve the least square problem min x r,im Ijx^ 1 ™ — A red x^™||. 

^r,im 

When the the matrix A rcd is not full column rank, we use the modified OMP 
algorithm which is adjusted to account for the errors due to noise. As noted 
above, in the noiseless case, one can expect a perfect reconstruction and thus 
the threshold error e can be set to or a very small number. However, with 
noisy signals, some care must be taken in choosing e. On the one hand, if the e 
is chosen too large, the algorithm may stop before a solution is reached. On the 
other hand, if e is chosen too small, the reconstructed signal may include bands 
that are not in the originating signal. The problem of too high threshold is 
solved by changing the stop criterion. Instead of stopping the algorithm when a 
threshold is attained, we check at each iteration whether the block that reduces 
the residual error the most causes the resulting matrix to be rank deficient. 
When this occurs, the iteration are stopped and the block is not added to the 
matrix. 

An additional change is made to the algorithm when treating the blocks. In 
the noiseless case each block corresponds to a single band in x[l]- When sampling- 
noisy signals, we divide each block into several sub-blocks. The reason for this 
division is that, with noisy signals, the identification of the bands is not accurate. 
Identified bands may be wider than the originating bands. This is particularly 
true if the threshold is chosen small. This widening may cause the inclusion of 
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Figure 6: Mean recovery times for equal 4-bands real signals for different spectral 
supports (total bandwidth) with i<Nyquist=40 GHz, total sampling rate Ftotai = 
12 GHz 
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Figure 7: Condition numbers histogram of the mixing matrix in (|21[) for 4 bands, 
200 MHz width real signals, F to taJ = 12 GHz 



^r,im 

false frequencies whose corresponding columns in A rcd are linearly dependent 
on the columns corresponding to the support of the originating signal. By using 
smaller sub-blocks such columns may be isolated from the rest of the columns 
in their corresponding band. 

The recovery scheme was tested against real-valued signals with 8 bands 
(4 positive frequencies bands and 4 negative frequencies bands). The signals 
without noise were generated and sampled exactly as in the noiseless simulations 
of real signals. Noise was added randomly at each frequency of the pre-sampled 
signal according to a normal distribution with standard deviation a = 0.04; 
the SNR was defined by 10 log 10 (l/((T-<jF max /F 2 ) ) = 10.5 dB. This definition 
takes into account the accumulation of noise in baseband due to sampling. The 
sampling rates were the same as those in the noiseless simulations. The indicator 
functions \ l [I] were constructed using the same parameters as those used in [9] . 
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Figure 8: Condition numbers mean value vs. number of aliased frequencies for 
4 bands, 200 MHz width real signals, F to tai = 12 GHz 
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Figure 9: Empirical success percentages for 4 bands of real signals noise that 
were contaminated with a noise with a standard deviation of <r = 0.04 for 
different spectral supports of the signal (-FLandau) with -Fkyquist^O GHz, and a 
total sampling rate Ftotai = 12 GHz 



Each band in x[l] was widened by 20 percent on each side. The sub-blocks used 
in the modified OMP had spectral width of 100 MHz. The success was measured 
by the algorithm's ability to achieve a low error li norm below 2a-J F max /F 2 = 
4.47<7 for each recovered band. The mean error for each recovered band X lec (f) 
and the true band X(f) were evaluated as follows: 

|4 J \X Iec (f) - *(f)\df < 2<T\/F~/F* 

where B is the band support. 

Statistics on recovering 8 bands 200 MHz width each are based on 10000 
tests. The simulation showed that, although the algorithm's performance in- 
evitably decreased, it still achieved a high empirical recovery rate (37 failures 
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out of 10000 tests). Additional simulations were performed by changing the 
Landau rates as was done in the simulations performed for the noiseless case. 
In Fig. [5] the empirical success percentage is presented for 1000 simulations of 
noisy signals. The results of the simulations are similar to those in the noiseless 
case. When the total sampling rate is 8 times higher than the Landau rate, 
high success percentage was achieved. The recovery error level depended on 
the threshold choice. Lower threshold allows more accurate reconstruction but 
increases the recovery time. Moreover, additional parameters adjustments are 
also necessary (widening percentage, sub-blocks size). Different error criteria 
are also possible. For example, choosing li norm instead of l\ norm and setting 
the error threshold to be 3.3a as in [S] resulted in 99.5 percent empirical success 
rate in recovering 1.6 GHz Landau rate signals and 99.8 percent for 1.5 GHz 
Landau rate signals. This is in contrast to simulations results in Fig. [9] where 
empirical perfect reconstruction was obtained for those signals under different 
criteria. 

5 Conclusions 

In this paper we describe a synchronized multirate sampling scheme for accu- 
rately reconstructing sparse multiband signals using a small number of sampling 
channels (3 in our simulations) whose total sampling rate is significantly lower 
than the Nyquist rate. 

Although the same data that is obtained from our scheme can be obtained 
by a multicoset scheme, such a multicoset scheme requires many more channels 
and a time accuracy that cannot be attained practically. Moreover, our scheme 
processes the data differently, in a way that results in significantly wider base- 
bands and thus greatly reduces the effects of aliasing. By synchronizing the 
sampling channels our scheme is able to reconstruct signals correctly even in 
most cases in which the signal is aliased in all the sampling channels. 

The main advantage of multicoset sampling schemes is theoretical. Because 
in a multicoset sampling scheme each channel samples at the same frequency, 
one has a mathematical structure that enables a perfect reconstructions of ideal 
multiband signals from samples taken at a total sampling rate that is close to a 
theoretical lower bound. This bound is attained only under special assumptions 
regarding the number and width of the signal bands. Moreover, the bound 
requires the number of sampling channels to be twice the number of signal 
bands |5j. Hence, in many cases the number of sampling channels becomes too 
high for a practical implementation. 

The main advantage of the SMRS is in the use of a small number of sampling 
channels that operate at relatively high rates that can reconstruct accurately 
sparse signal that consist of several bands. Sampling at higher rates has a 
fundamental advantage in that it increases the SNR after sampling. Another 
advantage is that implementation of the SMRS scheme does not require a priori 
knowledge of the maximum width of the signal bands. A third advantage of our 
scheme is that, in many cases the signal can be reconstructed by simple matrix 
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inversion rather than through a search algorithm as in the multicoset recovery 
scheme of [5]. 

In the SMRS scheme, when sampling sparse signals, most of the sampled 
spectrum is unaliased in at least one of the sampling channels. Hence, most 
of the spectrum can be reconstructed directly from the unaliased parts of the 
spectrum. On the other hand, in multicoset sampling schemes an alias in one 
channel is equivalent to an alias in all channels. Furthermore, a multicoset 
scheme downconverts signals to much lower frequencies, thus increasing the 
negative effects of aliasing. Our numerical simulations indicate that the reduced 
aliasing in our SMRS scheme results in significantly better performance over a 
multicoset sampling scheme of [5] whose number of sampling channels is small. 
Also, due to reduced aliasing, it is expected that the reconstruction in the SMRS 
scheme will be robust when noise is added in the sampling process. 

Although we do not have a rigorous criteria for perfect reconstruction, by 
examining the multicoset pattern that yields the same data it might be possible 
to obtain necessary conditions for a perfect reconstruction. 

Our scheme also solves ill-conditioned linear equations by using a modified 
OMP algorithm. Whereas we obtained satisfactory results with it, we do not 
have criteria for determining when the modified OMP algorithm converges to the 
correct solution. However, this shortcoming may be not as significant because 
of the possibility of using other algorithms to solve the equations. 

6 Appendix 

In this appendix we present the modifications to ([5]) for the real signals recovery 
Since the signal is real-valued, its spectrum fulfills 

X(f)=X(-f) (14) 

where a + bj = a — bj is the complex conjugate and a and b are real numbers. 

It follows from (fl"4f and (fTJ), that for each channel index i, all the information 
about X l (f) is contained in the interval [0, F l /2]. Consequently, it is convenient 
to choose the sampling frequencies F l such that F % j2 = A/M l /2 where M l /2 
is an integer. Because the conjugation operation a + jb : a + jb i— > a — jb is not 
complex linear, §5§ needs to be replaced with two systems of equations; one for 
the real part and one for the imaginary part. 

We use the following notations to represent the spectrum of the real signals 
in the discretized frequencies: 

X l [k] =X\kAf) k = -[M i /2j,...,[M i /2\, (15) 
X[k] = X(kAf) k = - [M/2\ , . . . , [M/2J . 

The sequence X l [k] contains the samples oiX z (f) in the baseband [— F % /2, F l /2\. 
The sequence X[k] contains the samples of X(f) given in [— MA//2, MA//2], 
where M is chosen to fulfill M = [.Fkyquist / A/] . Equation ([3]) now takes the 
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following form: 

LM/2J oo 

X l [k]=F l X W J2 {k + nM r )}. (16) 

l=—\M/2\ n=~ca 

Equation (|16p can be written in a matrix form as 

x 1 = A»'x (17) 

where x 1 and x are given by 

(x l ) fe+L AfV2j+i = ^[k], - [M i /2\ < k < LMV2J, (18) 
(x) fc+L M/2j+i =X[k], - LA//2J < k < [M/2\, 

and A* is a matrix whose elements are given by 

oo 

AU L M./2j+ u+ LM/2j + i=^ J2 5[l-{k + nM% (19) 

n— — oo 

Note that, since the signal is real valued, all of its spectral information is con- 
tained in the positive frequencies. 

Each element in A 1 is equal to either F l or 0. Equation (jTTJ) for the different 
sampling channels can be concatenated as in complex signals case to yield 

x = Ax. (20) 

The spectrum can be decomposed into its real and imaginary parts. As a result 
(|2"D|) becomes 

x r = A" x r , (21) 
x™ = A "V™ 



where x r = Re(x) and x lm = Im(x). In addition only components that cor- 
respond to positive frequencies are retained in the vectors x r and x l . The 
elements of the matrices A and A are given by 

Al, LM /2j_i+i = A M+1 + A fc|M _,, 1 = 0,..., [M/2\, (22) 
A fe,LAf/2j-;+i = A fe,M-; - Afc,j + i, 1 = 0,..., [M/2\. 
The reconstruction is performed with (f2"Tj) exactly as in the complex case. 
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